--- title: Vertex Dynamics Models keywords: fastai sidebar: home_sidebar summary: "Simulation codes" description: "Simulation codes" ---
device=torch.device('cpu')
if torch.cuda.device_count():
device=torch.device('cuda:1')
dtype = torch.float32
import networkx as nx
# let's seed RNG for sanity and reproducibility
np.random.seed(42)
# define cell monolayer
v_x,regions =unit_hexagons(10,10) # unit hexagons
# convert Voronoi regions to cells and edges
edge_list,cells = VoronoiRegions2Edges(regions)
# perturb vertices
v_x += np.random.randn(v_x.shape[0], v_x.shape[1])*.2
cell_graph = Monolayer(vertices=Vertex(v_x.copy().tolist(),dtype=dtype),
edges=torch.tensor(edge_list),
cells=cells)
cell_graph.to_(device)
Gnx,pos=graph2networkx_with_pos(cell_graph)
fig = plt.figure(figsize=[5,5])
fig.clf()
ax = fig.subplots()
ax.axis(False);
nx.draw(Gnx,pos,node_size=50,width=4,ax=ax,node_color='#FF00FF',edge_color='#51C5FF')
plt.show()
plt.close()
# Define energy function
omega = torch.tensor([np.pi/2],dtype=dtype,device=device)
phase = torch.rand((len(edge_list),1)).type(dtype).to(device)
p_keep = 0.2 # fraction to keep active
e_ij_on = (torch.rand((len(edge_list),1))<p_keep).type(dtype).to(device)
def monolayer_energy(monolayer_,Tau):
Perm = cell_graph.perimeter()
Area = cell_graph.area()
Leng = cell_graph.length()
direction = cell_graph.direction()
dir_coeff = torch.abs(direction[:,1])/torch.norm(direction,dim=1)/5
gamma_ij = e_ij_on*dir_coeff.view(-1,1)*torch.cos(omega*Tau+phase)**2
gamma_ij_by_lij = Leng*gamma_ij
Energy_tot = torch.sum(.01*(Perm)**2) + torch.sum((Area-2.3)**2) + torch.sum(gamma_ij_by_lij)
return Energy_tot, gamma_ij
# Networkx's edge ordering is different
edge_idx = dict(zip([tuple(e) for e in cell_graph.edges.tolist()],range(cell_graph.edges.shape[0])))
edge_idx_order = [edge_idx[e if e in edge_idx else (e[1],e[0])] for e in Gnx.edges ]
def draw_w_tension(i):
pos = dict(zip(range(verts_t[i].shape[0]),verts_t[i].numpy()))
ax.cla()
ax.axis('off')
ax.set_title(f'Time : {verts_frames[i]:2.3f}')
#node_color=range(24), node_size=800, cmap=plt.cm.Blues
nx.draw(Gnx,pos,node_size=50,width=4,ax=ax,node_color='#FF00FF',
edge_color=line_tensions[i].numpy().squeeze()[edge_idx_order],edge_cmap=plt.cm.bwr)
_,Contractions = monolayer_energy(cell_graph,0)
fig = plt.figure(figsize=[7,7])
fig.clf()
ax = fig.subplots()
ax.axis(False);
nx.draw(Gnx,pos,node_size=50,width=4,node_color='#FF00FF',
edge_color=Contractions.cpu().numpy().squeeze()[edge_idx_order],edge_cmap=plt.cm.bwr)
# Simulation param-s
Dt = 2**-5 # time step size
t = [0]
Energies = []
Forces = []
verts_t =[]
verts_frames=[]
line_tensions=[]
t_total = 2**10
cell_graph.vertices.requires_grad_(True)
# Numerical integration
print('Integration (Euler\'s method):')
for n in range(t_total):
cell_graph.set_zero_grad_() # reset grad accumulator
t.append(t[-1]+Dt) # update last frame time
# total potential energy of the system:
E,Contractions = monolayer_energy(cell_graph,t[-1])
Energies.append(E.item()) # E(t-1)
E.backward() # compute gradients
dxdt = -cell_graph.get_vertex_grad() # dx/dt=-dE/dx
# Update vertex position
with torch.no_grad():
cell_graph.vertices.x += dxdt*Dt
Forces.append(torch.norm(dxdt,dim=1).mean().item())
verts_t.append(cell_graph.vertices.x.detach().cpu().clone())
verts_frames.append(t[-1])
line_tensions.append(Contractions.detach().cpu().clone())
if round((n+1)%(t_total/8))==0:
mean_grad = torch.norm(dxdt,dim=1).mean().item()
print(f't={t[-1]:2.3f}: E={E.item():5.4g}; aver |dx/dt|= {mean_grad:3.2g}')
Energies.append( monolayer_energy( cell_graph,t[-1])[0].item() )
plt.figure(figsize=[5,3])
plt.plot(t,Energies);plt.xlabel('time');plt.ylabel('energy');
# add forces (except last frame)
ax2=plt.gca().twinx()
ax2.set_ylabel('Average $|F|$',color='red')
ax2.plot(t[:-1],Forces,'r-',alpha=.4,lw=2)
plt.show()
# Print final Perimeters and Areas
print(f"Perimeters:{cell_graph.perimeter().detach().squeeze()}\nAreas:{cell_graph.area().detach()}")
import matplotlib.animation as animation
from IPython.display import HTML
import PIL
fig = plt.figure(figsize=[7,7])
fig.clf()
ax = fig.subplots()
ax.axis(False);
draw_w_tension(0)
plt.show()
plt.close()
anim = animation.FuncAnimation(fig, draw_w_tension, interval=200,
frames = range(0,len(verts_t),max(1,round(len(verts_t)/64))))
HTML(anim.to_jshtml())
Save as GIF
# len(verts_t)
pil_images=[]
for i in range(0,len(verts_t),5):
fig = plt.figure(figsize=[7,7])
fig.clf()
ax = fig.subplots()
ax.axis(False);
draw_w_tension(i)
# plt.show()
canvas = plt.get_current_fig_manager().canvas
canvas.draw()
pil_images.append(PIL.Image.frombytes('RGB', canvas.get_width_height(),
canvas.tostring_rgb()) )
plt.close()
pil_images[0].save('sim_anisotropic_contraction.gif', format='GIF',
append_images=pil_images[1:], save_all=True, duration=100, loop=0)